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Abstract 

An efficient semi-implicit second-order-accurate finite-difference metfiod is 
described for studying incompressible Rayleigh-Benard convection in a box, 
with sidewalls that are periodic, thermally insulated, or thermally conducting. 
Operator-splitting and a projection method reduce the algorithm at each time 
step to the solution of four Helmholtz equations and one Poisson equation, 
and these are are solved by fast direct methods. The method is numeri- 
cally stable even though all field values are placed on a single non-staggered 
mesh commensurate with the boundaries. The efficiency and accuracy of the 
method are characterized for several representative convection problems. 
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I. INTRODUCTION 



Experiments over the last three decades have discovered many fascinating and poorly 
understood examples of pattern formation in large-aspect-ratio Rayleigh-Benard convec- 
tion [1,2]. Because of the prominent role that these experiments play in understanding 
sustained nonequilibrium systems [2] and because many of the observed phenomena such 
as spatiotemporal chaos are difficult to analyze mathematically [3], there is a need to de- 
velop computer codes that can simulate these experiments quantitatively so that theory and 
experiment can be compared with one another. Once validated, such codes can further be 
used to explore regimes not easily attained by experiment such as low Prandtl number, and 
to calculate quantities that are difficult to deduce from experimental data, such as mean 
flows [4,5] and fractal dimensions [6]. 

The regime of large aspect ratio F (ratio of horizontal fluid width to fluid depth) poses 
significant computational challenges. Many numerical degrees of freedom (basis functions or 
mesh points) are needed to represent the spatial features of the fluid and often the dynamics 
needs to be studied over long times (many multiples of the horizontal thermal diffusion 
time Th = r^r^, where = (f/n is the vertical thermal diffusion time defined in terms 
of the fluid depth d and fluid thermal diffusivity k) to insure that nontransient behavior 
is being observed. Since the largest time step allowed by numerical stability for exphcit 
or semi-implicit algorithms (the ones most commonly used in Navier-Stokes calculations) 
is typically 0.05 Ty or smaller, simulations in a representative F = 50 cell may require 
10^ or more time steps to eliminate a transient and then study the statistically stationary 
properties of the asymptotic dynamics. The many degrees of freedom and long integration 
times together with the need to repeat runs for different parameter values imply that efficient 
algorithms are essential for studying the large-aspect-ratio regime. 

Because of these computational challenges, there have been few simulations of three- 
dimensional Rayleigh-Benard convection with aspect ratios exceeding 10. Recent calcula- 
tions with F as large as 64 have been carried out by Pesch and collaborators, who used a 
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pseudospectral code on serial and parallel computers to study spiral defect chaos, rotating 
convection, and other problems [6-8]. However, their code uses periodic boundaries and so 
can not take into account quantitatively the influence of lateral walls on the bulk dynamics. 
Arter and Newell [9] and Tomita and Abe [10] have carried out simulations in large boxes 
with thermally insulated no-slip sidewalls, the former in a 16 x 11.5 aspect ratio box, the 
latter in a F = 18.84 square box. Xi et al [11] have studied the transition to spatiotemporal 
chaos of a convecting fluid in a F = 60 square cell, but with free-slip horizontal boundaries 
that are difficult to achieve experimentally. Finally, a Caltech-Duke collaboration has re- 
cently reported results [12, 13] obtained with a parallel spectral element code [14] for aspect 
ratios up to 30. Their code can treat quantitatively most geometries and lateral boundaries 
used by experimentalists, including ramps [15], spoiler fins [16], and lateral walls of finite 
thickness and finite thermal conductivity. However, the generality of the spectral element 
algorithm makes it substantially more expensive to run than algorithms optimized for a 
simple geometry such as a box or cylinder. 

In this paper, we introduce and analyze an efficient semi-implicit finite-difference algo- 
rithm for studying incompressible Rayleigh-Benard convection in a box, with lateral walls 
that are periodic, thermally insulated, or thermally conducting. The code complements 
the more fiexible spectral element approach [12, 13] by being more than an order of mag- 
nitude more efficient on a serial processor, for a box with these boundary conditions. It 
is well suited for studying long-time dynamics of small- to moderate-aspect-ratio boxes [4] 
(F < 20), with lateral boundaries that are close to those of many experiments, although not 
fully quantitatively accurate since the finite thickness and finite thermal diffusivity of the 
lateral walls is not taken into account. 

The main advantages of our algorithm are the simplicity of implementation and its 
efficiency on a single processor. Its simplicity arises from the use of a single mesh for all field 
values. (This is called a "non-staggered" or "collocated" mesh in contrast to a "staggered" 
mesh for which the values of different fields appear at different points in space [17-19]). A 
non-staggered mesh reduces the effort to write and to validate a code (compared to one 
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using a staggered mesh), and facilitates porting the code to a distributed- memory parallel 
computer. Earlier work on Navier-Stokes integrators has suggested that non-staggered mesh 
codes can be numerically unstable because of pressure oscillations [17, 20]. Our results below 
show that an algorithm to integrate the Boussinesq equations on a non-staggered mesh can 
be numerically stable. 

The use of a single non-staggered mesh also helps to explain the efficiency of the algo- 
rithm. Using a standard operator splitting and projection method together with second- 
order-accurate finite differences [21,22] on a uniform three-dimensional mesh, the advance- 
ment of the velocity, temperature, and pressure fields at each time step requires the numer- 
ical solution of four Helmholtz equations and one Poisson equation. Because these elliptic 
equations and their boundary conditions are separable, they can be solved efficiently using 
fast direct methods from the FISHPACK library [23,24], with a complexity per problem 
of 0(A'"log(A^)), where N is the total number of mesh points. Fast direct methods are more 
efficient than most iterative methods on a single processor [25], and have the additional 
advantage that no internal parameters need to be adjusted to obtain convergence. However, 
fast direct methods are not applicable to complex geometries, to problems with spatially 
varying parameters, or to complicated boundary conditions that lead to nonseparable equa- 
tions. 

The remainder of this paper is organized as follows. In Section II, we discuss details of 
our algorithm, namely how the fields and equations are discretized and how the resulting 
equations are solved. In Section III, we discuss the convergence properties of the algorithm 
and its efficiency for several representative two- and three-dimensional convection problems. 
We confirm empirically the second-order accuracy of the solution and examine how the 
largest time step allowed by stability varies with Prandtl number and with Rayleigh num- 
ber. Finally, Section IV presents our conclusions and suggests some avenues for further 
algorithmic improvements. Applications of the algorithm to study quasiperiodic dynamics 
and spiral defect chaos in three-dimensional boxes can be found in Refs. [4, 5] . 
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II. DETAILS OF THE ALGORITHM 



A. Equations and Boundary Conditions 

Our goal is to integrate the Boussinesq equations that describe incompressible buoyancy- 
driven Rayleigh-Benard convection with an external force f . These equations can be written 
in the dimensionless form [26] 

dtT{t,x,y,z)^[-{v'V)T] +V'T, (1) 
dtv{t, X, y, z) = [-(v • V)v + aRTz + erf] + crV^v - Vp, (2) 
V'v = 0. (3) 

The variables x and y denote the horizontal coordinates, while the z variable denotes the 
vertical coordinate, with the unit vector z pointing in the direction opposite to the grav- 
itational acceleration. The field v = {vj:,Vy,Vz) is the velocity field at point {x,y,z) at 
time t, while p and T are the pressure and temperature fields respectively. The dimen- 
sionless parameters a and it! denote the Prandtl and Rayleigh numbers respectively. The 
vector field f (i, x, T, v) is some external force, e.g., a Coriolis force f 2v x f2 arising from 
a rigid rotation of the convection cell with constant angular velocity CI = Qz. The terms 
grouped in brackets in Eqs. (1) and (2) are those containing nonlinear terms or linear terms 
with low-order spatial derivatives, and will be integrated explicitly by the operator splitting 
method described below. 

We would like to integrate Eqs. (l)-(3) in a box geometry defined by the region 

2 - - 2 ' 2 - 2 ' 2 - - 2' ^ ^ 

where and Ty are the aspect ratios in the x and y directions respectively (the depth of 
the fluid has length 1). A no-slip velocity condition on all material walls is assumed 

V = 0, (5) 

and the temperature T is constant on the bottom and top plates. 




for 



z 




1 



(6) 



The code allows the temperature on any opposing pair of lateral walls to be periodic, or 
to satisfy on each lateral wall an arbitrary Dirichlet boundary condition (e.g., a thermally 
conducting wall corresponding to a linear conducting profile of the form T — a + bz, where a 
and b are constants), or an arbitrary Neumann condition (e.g., a thermally insulating wall 
with dnT = 0, where dn is the normal derivative to the boundary at a given point). To 
simplify the following discussion, we will consider only the case of insulating sidewalls 



since the other cases involve just simple modifications. Although the pressure field p formally 
has no associated boundary condition since it does not satisfy a dynamic equation, we will 
be imposing a Neumann boundary condition on p as we explain below (see Eq. (17)). 



We next discuss the time integration method, since its structure can be explained before 
having to specify a spatial representation for the fields. In the following subsection, we 
discuss how the fields and equations are discretized and the latter solved using second- 
order-accurate finite differences on a uniform spatial mesh. 

Our time integration method uses a standard operator splitting and projection 
method [21, 22], in which the nonlinear terms containing lower-order or no spatial derivatives 
are integrated explicitly, then the linear diffusion operators are integrated implicitly, and 
finally the pressure term — Vp is integrated to project the velocity field at the next time step 
into the space of divergence-free velocity fields. Operator sphtting has two benefits. First, 
the evolution equations for T and for the velocity components are decoupled from one 
another, which simplifies the overall algorithm and substantially reduces the total computer 
memory needed. Second, operator splitting allows larger time steps since the largest time 



dnT = 0, 



on lateral walls. 



(7) 



B. The Time Integration Method 
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step At allowed by stability is bounded by a first power of the spatial resolution Ax, rather 
than by a second power as would be the case for a fully explicit method. 

Let us assume that, at the nth time step t„ = nAt with n > 0, initial fields T"' and v" 
are known that are consistent with the boundary conditions Eqs. (5)-(7). These fields are 
then advanced to the future values T"+^ and v"+^ at time tn+i — tn + At as follows: 

1. The nonlinear advective term Nt[T,v] — — (v V)T of Eq. (1) is integrated explicitly 
using a second-order-accurate Adams-Bashforth method 

T* = T'^ + ^(3iVT[T",v"]-iVr[T"-\v"-^]), n > 0, (8) 

to produce an intermediate field T*. Here T"~^ and v"~^ denote field values stored 
from the previous time step tn-i — tn — At. For the first time step n — only, a second- 
order-accurate single-step integrator, Heun's method [27], is used in place of Adams- 
Bashforth to avoid the dependence on the unavailable field values at time t — —At. 

2. The intermediate field T* is then advanced to the temperature field T"+^ at time tn+i 
by using as initial data for an imphcit Crank- Nicolson step apphed to the diffusion 
term in Eq. (1): 



At 



- (v2r"+^ + V^T") . (9) 



This can be written as a constant-coefficient Helmholtz equation for the future 
field r"+^: 

_ ^v^) r"+i = T* + ^V^T", (10) 

and is solved with the boundary conditions Eqs. (6) and (7) applied to T"+i. 

3. Three similar pairs of explicit and implicit steps are then executed successively, first 
for the velocity component v.j., then for Vy, and then for v^- The explicit steps advance 
the velocity field v** to an intermediate field v*, and then the implicit steps advance v* 



to a second intermediate field v**. If we denote by Ni[T, v] the expressions in brackets 
of Eq. (2) for i — x, y, and z, then each exphcit step has the form 

< = < + ^ (3^^ [T, v1 - N, [t, v"-i] ) . (11) 

Heun's single-step method is again used at time to = to avoid unavailable field values 
at time t_i = —At. Each field f " is next used as initial data for an implicit Crank- 
Nicolson step that yields a constant-coefficient Helmholtz equation for the field v**: 

1 _ -^v') vr = < + ^v^„r. (12) 

This is solved with the no-slip boundary condition v** = on all surfaces, Eq. (5). 

4. An incompressible velocity field v"+^ at time t^+i is obtained from the field v** by 
integrating the final operator step 

dtv = -Vp, (13) 

with initial data v**, followed by a projection method [22,28]. We approximate the 
time derivative in Eq. (13) with a first-order- accurate stencil. 



v"+i - V** 



At 



_1 (Vp-+i + Vp"), (14) 



apply the divergence operator to both sides, and then use Eq. (3) in the form 
y ,yn+i _ Q rpj^-g yjg^g Poisson equation for the pressure field p: 

V V^' = -V V + ^ V • V**. (15) 

Once p is known by solving Eq. (15), we obtain v"+^ from Eq. (14) in the form 

At 



(Vp"+^ + Vp"). (16) 



Although mathematically there is no boundary condition for p — and by discretizing 
first space and then time, a boundary condition for p can be avoided as explained in 
Refs. [21,29,30] — we will solve Eq. (15) with the Neumann condition 
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= on all walls, 



(17) 



since this allows us to use a fast direct method to solve Eq. (15). There is a substan- 
tial literature concerning the appropriateness and accuracy of the boundary condition 
Eq. (17) [31, 32]. Rather than review this literature, we simply point out that a Neu- 
mann pressure boundary condition has been shown by previous researchers to produce 
acceptably accurate results for problems in which the fluid is confined by no-slip sur- 
faces, and we show directly in Section III that our algorithm is second-order accurate 
in space and second-order accurate in time for several representative problems. 

The most time consuming part of this algorithm is, by far, solving the four Helmholtz 
equations Eq. (10) and Eq. (12) (for i = x ,y, and z) and solving the Poisson equation 
Eq. (15). 



C. Discretization on a Uniform Mesh 

The explicit and implicit steps of the previous section — Eqs. (8) and (10), Eqs. (11) 
and (12), and Eqs. (15) and (16) — are carried out by discretizing the fields and equations 
on a single non-staggered mesh of points 

Xyfc = (iAxjAy, kAz) , (18) 

that is commensurate with the sides of the box Eq. (4). The mesh indices (i, j, k) satisfy 

-^<i<^ _^<o<^ -^<k<^. (19) 
2--2' 2~''^~2 2~~2 ^ ^ 

The aspect ratios (F^;, Ty) and the positive integers {N^, Ny, N^) are specified as input to the 
code, and the corresponding spatial resolutions (Ax, Ay, Az) are then determined from the 
relations Ax — V^/N^, Ay = Ty/Ny, and Az = l/{Nz). For large-aspect-ratio convection 
problems, typically Ax — Ay > Az since the x and y directions are equivalent and there is 
a finer structure in the vertical direction caused by the close opposing horizontal plates. 
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At all mesh points Eq. (18) interior to the box, the first- and second-order spatial deriva- 
tives are approximated using centered second-order-accurate 3-point finite-difference sten- 
cils. If (xjjfc) denotes the values of a field m(x) at the mesh points, then the partial 
derivative d^u at Xj^fe is approximated by 

with similar expressions for dyU and d^u. The Laplacian of u at Xy^ is approximated by the 
usual 7-point stencil 



ijk ~ Ax^ Ay2 ^ 

+ Az' ■ ^^^^ 

Nonsymmetric finite-differences are needed to evaluate expressions on those boundaries 
for which a Neumann condition holds (we will call these "Neumann boundaries" ) since field 
values outside the domain are not available. Thus the right side of the Helmholtz equation, 
Eq. (10), needs to be evaluated on the Neumann boundaries for which Eq. (7) holds. The 
value of V^T* can be approximated there to second-order accuracy by using one-sided 4- 
point finite-difference approximations for the 2nd-order derivatives, e.g., 

v^x^ \ojk ' y^^) 

with similar expressions for d^T* at x = Vx, and for dyT* on the y = ^ and y = Vy 
boundaries. The divergence V»v** on the right side of the pressure equation, Eq. (15), can 
be approximated to second-order accuracy using the Dirichlet data Eq. (5) and interior field 
values by replacing Eq. (20) with the following 3-point one-sided finite difference 

r;^ 1 _ -3(^x)0jfc + '^{Vx)ljk - {Vx)2jk 4^{Vx)ljk ' {Vx)2jk 

[dxVxhk ^ ^ , (24) 

with similar expressions for dyVy and d^Vz- The advective derivative —{v'V)T in Eq. (8) 
vanishes on these Neumann walls since v does, and so the explicit steps do not require 
special treatment. 
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Given the discretizations Eqs. (20)-(24), the exphcit time steps Eqs. (8) and (11) are 
easily evaluated at all interior points and on the Neumann boundaries. For the implicit 
steps, the right sides of Eqs. (10) and (12) are also evaluated on the interior mesh points and 
on the Neumann boundaries. These right sides are then used as input to the FISHPACK [24] 
fast direct solvers hwScrt in three dimensions or hwscrt in two dimensions. Also provided 
as input to the FISHPACK solvers are the corresponding boundary conditions, Eqs. (6) 
and (7) for T, Eq. (5) for the velocity components, and Eq. (17) for p. The FISHPACK 
solvers return second-order- accurate values (with respect to the spatial resolution) of T, v, 
and p on the mesh Xjjfc. 

We conclude this section with the observation that the discrete velocity field v""*"^ ob- 
tained from the concluding step Eq. (16) is only approximately divergence-free even on the 
mesh points Xj^fc, i.e., V • v**"*"^ = 0{h?) where h is the larger of the spatial resolutions Aa;, 
Ay, and Az. This is because the discrete approximation Eq. (20) for the pressure gradient 
in Eq. (16) is not consistent with the discretization Eq. (21) used to approximate the Lapla- 
cian V* V on the left side of Eq. (15). The discrete Laplacian can be considered as arising 
from the evaluation of pressure gradients from pairs of nearest neighbor points as follows 



In contrast, Eq. (16) evaluates the pressure gradient at a point using a finite difference 
Eq. (20) that spans three mesh points. 



In this section, we discuss several tests that quantify the accuracy of the above algorithm 
for a convecting fluid in a two-dimensional rectangular domain with periodic sidewalls and 
in a three-dimensional rectangular domain with perfectly insulating sidewalls, Eq. (7). We 
first confirm the second-order accuracy of the code with respect to the spatial and time 




(25) 



(26) 



III. ACCURACY AND EFFICIENCY OF THE ALGORITHM 
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resolutions by studying how the temperature and velocity fields converge with increasing 
spatial and time resolutions respectively. We next show empirically how the maximum 
stable time step varies with the Rayleigh number R and the Prandtl numbers a. We then 
calculate the critical Rayleigh number Rc and plot the Nusselt number N{R) as a function 
of the Rayleigh number R, and obtain good agreement with an analytical expression [34] 
and with a spectral code [35]. Finally, we show the spatial structure of the fields near onset, 
to allow comparison with experiment [36] and with other codes. 

We note that, on a workstation with a 667 MHz 21264A 64-bit Alpha processor, a square 
box with aspect ratio F = 40 and spatial resolution Ax = Ay = Az = 1/8 takes about 4.8 s 
per time step of At — 0.001 1^. This corresponds to 80 minutes per vertical diffusion time ty 
and 90 days per horizontal diffusion time th, so this code is too slow to explore F > 20 cells 
over time scales exceeding a horizontal diffusion time. We discuss two ways of improving 
the efficiency of the code in our concluding comments of Section IV. 

A. Second-Order Convergence With Respect to the Spatial and Time Resolutions 

We begin by showing that the order of convergence p of the code is asymptotically 2 
(second order) in the limits of sufficiently fine spatial and time resolutions. By definition, 
the convergence with respect to spatial resolution is of order p ii \\uh — itexactH = 0{h^) in 



the limit /i 0, where ||-u|| = denotes the Euclidean norm of a field u on the 

spatial mesh, h — Ax — Az is the uniform spatial resolution in the x and z directions of a 
two-dimensional box, Uh{x, z) denotes a discrete numerical field on a mesh of resolution /i, 
and Ucxa_ct{x^z) is the unknown exact field on the spatial mesh. By writing Uh{x,z) = 
Inexact (2^, z) + C{x, z)h^ in the limit /i ^ 0, for some function C independent of /i, we deduce 
that the order p can be estimated by examining the quantity 



in the limit /i 0. The estimate Eq. (27) involves field values at the three levels of 
resolution 4/i, 2/i, and /i, coarsest to finest. A similar definition for the order of convergence 





(27) 
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with respect to time resolution can be made if the spatial mesh h is replaced by the time 
step g = At, 



We first studied the convergence with respect to the spatial mesh h for a two-dimensional 
box with periodic sidewalls, for parameter values — Stt/^c = 2.016, it! = 1725.0 1.01it!c, 
and a — 0.71. The initial conditions consisted of a small random perturbation about the 
linearly conducting state To = —z, vq = {uq, wq) = 0, and these were integrated until a 
stationary state was attained consisting of two rolls at the critical wave number Qc- For this 
small cell, an integration time of St^ was sufficient for the dynamics to become stationary. 
We then studied the temperature field, Th{x, z), and the x-component of the velocity field, 
Ufi{x, z), for different spatial resolutions N — T^/h — 16, 32, 64, and 128. The time step At 
was set respectively to the values At = 0.01, 0.005, 0.0025, and 0.00125 since the operating 
splitting makes the largest stable time step proportional to h. Table I summarizes the values 
of the hmit Eq. (27) and shows that indeed p/i — > 2 as /i — > 0, i.e., the code is asymptotically 
second-order accurate with respect to the spatial resolution h. 

We have also studied the convergence with respect to the time step g for a three- 
dimensional box with perfectly insulating sidewalls and for parameter values — Ty — 2, 
R — 1725.0 ~ 1.01i?c, and a = 0.71. The initial condition consisted of small random 
thermal perturbations, and these were integrated up to 20 diffusion times at which point 
the state became stationary. For various time resolutions g — At — 0.0001, 0.00005, and 
0.000025, all with a space resolution oi N — 64, the convergence was found (using Eq. (28)) 
to be p = 1.68. This provides evidence that the code is indeed asymptotically second-order 
accurate with respect to the time resolution g. 

B. Dependence of Maximum Stable Time Step on Rayleigh and Prandtl Numbers 

Since an important practical feature of any production code is the largest time step 
that can be taken before numerical instability occurs, we have studied the maximum stable 




(28) 
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time step as a function of the Rayleigh and Prandtl numbers. A three-dimensional box with 
periodic sidewalls and aspect ratio — Fy — 2 was used, with a spatial resolution Ax = 16. 
The Euclidean norm of the temperature field, ||T||, was calculated for various values of At 
each time after a interval of 20 vertical diffusion times so that transients decayed. The 
maximum stable time step was then defined as the value of At such that ||T|| remains 
bounded, i.e., ||T|| < 10^. 

In Fig. 1(a), we plot the maximum stable time step as a function of Rayleigh number 
for the two Prandtl number values a = 1 (square symbols) and cr = 10 (cross symbols). We 
see that the maximum stable time step decreases rapidly with increasing Rayleigh number. 
This is to be expected since the magnitude of the velocity and temperature fields increase 
with increasing R. In fact, a best log-log fit to the data yields the relation 

max(Ai) oc i?" (29) 

where a — —1.2 when a — 1, and a — —1.3 when a — 10. 

In Fig. 1(b), we plot the maximum stable time step as a function of Prandtl number 
for fixed Rayleigh numbers it! = 2048 (square symbols) and R — 8192 (cross symbols). 
For R — 2048, the maximum stable time step decreases toward both small and large Prandtl 
numbers. For R = 8192, the maximum stable time step decreases toward small Prandtl 
numbers but is approximately constant at large Prandtl numbers. The smaller time step 
needed at small Prandtl numbers can be attributed to the more dynamical nature of the 
convective flow at small Prandtl numbers, such as the presence of spiral defect chaos [33] . 

C. Estimate of the Critical Rayleigh Number Rc 

A linear stability analysis of the Boussinesq equations about the linearly conducting 
proflle between two infinite horizontal no-slip plates shows that the critical Rayleigh number 
Rc ~ 1707.76 with critical wave number Qc ~ 3.117, and that the values of Rc and Qc are 
independent of the Prandtl number a [2] . We tested these predictions and so validated the 
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code by using a two-dimensional box of aspect ratio — '^tt/qc — 2.016, with periodic 
sidewalls, for Prandtl number a — 0.71. We used a uniform spatial resolution h — 1/N — 
Ax = Az and varied the number N of mesh points. 

The critical Rayleigh number Rc was estimated as the approximate value of R for which 
the growth rate A = A(i?) of a small- amplitude (0.01) random perturbation about the linear 
profile interpolated to zero as a function of R. Thus for a sufficiently tiny initial perturbation 
of the conducting profile, there is a time interval over which the velocity component w 
grows approximately exponentially 



where A is the growth rate, and c is independent of t but can vary with R. For R> Rc, the 
growth rate is positive, for R < Rc, the growth rate is negative and interpolating between 
known positive and negative values provides an estimate of Rc, for which A = 0. 

Our protocol was to set R — i?+ = 1725 — 1.01 i?c just above onset, set the initial 
velocity field to zero, vq = {uo,Wq) = 0, and set the initial temperature field Tq = —z + 6T 
to a tiny random perturbation ST{x, z) of the hnear profile T — —z, with \ST\ < 0.01. The 
initial conditions were then integrated for a short time and the growth rate estimated from 
the formula 



where t2 and ti < ^2 ai'c two times during the exponential growth of the magnitude of 
the z-component of the velocity field w. The calculation was then repeated with the same 
initial condition but for R — R_ — 1691 = 0.99i?c to estimate a decay rate A_. The critical 
Rayleigh number was then estimated as the zero of the line joining the points {R+, A+) 
and (i?_,A_). The estimated critical Rayleigh numbers Rc as a function of the number 
of mesh points N are summarized in Table II. The values are correct to a relative error 
of better than one percent for the finest spatial resolution, confirming the correctness and 
convergence of the discretization and of the solution technique. 




(30) 



A 



lR{\\w{t2,X,z)\\/\\w{t,,X,z)\\) 
t2 — ti 



(31) 
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D. The Nusselt number versus Rayleigh number Curve N{R) 



Another way to characterize the accuracy of a convection code is by the dimensionless 
Nusselt number N{t, R, a), which is the instantaneous global vertical heat transport through 
the fluid layer, normalized to the heat transport arising from thermal conduction alone. For 
the dimensionless variables used in Eqs. (l)-(3) above, the Nusselt number can be expressed 
in the form [26] 

7V=l + (^(T-Teond)), (32) 

where w is the z-component of the velocity field and Tcond = —z is the temperature profile 
of the linear conducting state with v = 0. The brackets (• • •) denote an average of a 
quantity over the horizontal coordinates. Sufficiently close to onset, numerical values of N 
can be compared with an analytical expression [34] that is valid asymptotically in the limit 
R-R^^ 0+. 

We have evaluated Eq. (32) for the two-dimensional domain with periodic sidewalls of 
the previous section, with parameters = 2.016, a = 0.71, and N = 16. Starting from a 
small perturbation of the linear profile, we integrated until a stationary state was attained, 
and then evaluated Eq. (32) for the stationary state. Fig. 2 shows how N empirically 
varies with R, and we compare this curve with the analytical result of Schliiter et al [34], 
and with numerical values obtained by Clever and Busse[35], who used a two-dimensional 
spectral code with periodic side walls. The agreement is good in both cases and confirms 
the correctness and accuracy of the code. 

E. Spatial Structure of the Numerical Solutions 

We conclude this section with a few examples of the spatial structure of the fields obtained 
from the code, to show that the fields are physically reasonable when adequately resolved, 
and are qualitatively in agreement with other codes and with experiment [36]. 
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for (7 = 0.71, Fig. 3 shows contours of constant temperature T and for two values of the 
Rayleigh number, R — 2500 in (a) and R — 10^ in (b). Warm fluid ascends in the middle 
of the cell and descends as cooler fluid on both sides. As R increases, a thermal boundary 
layer forms at the top and bottom plates, creating a finer spatial structure that will require 
eventually a decrease in the vertical spatial mesh size iS.z. Fig. 4(a) shows the corresponding 
velocity field v = (m, w), while Fig. 4(b) shows the vertical component w through the midline 
of the cell. The occurrence of two square-shaped convection cells of opposite vorticity is in 
good agreement with experiment [36]. 

Fig. 5 shows constant temperature contours in a three-dimensional box with insulating 
sidewalls at time t — 200 i^, for parameters F = 16, i? = 2500, a — 0.71, and h — Ax = 
A?/ — Az = 1/8, At = 0.01. In agreement with experiment [38] and with calculations on the 
Swift-Hohenberg model of convection [39], the rolls are approximately normal to the lateral 
walls and the pattern consists of two diagonally opposite foci. For shghtly higher R — 
8500, Fig. 5(b) shows that the oscillatory instabihty commenced in the form of ripples 
that propagate along the length of the rolls. The occurrence of the oscillatory instability 
and its spatial form are in good agreement with the linear stability analysis of Busse and 
collaborators [26] and with experiment. 

IV. CONCLUSIONS 

We have described and characterized a semi-implicit finite-difference algorithm for inte- 
grating the Boussinesq equations in two- and three-dimensional boxes, with sidewalls that 
are periodic, thermally insulated, or thermally conducting. Our approach is useful for sim- 
ple geometries like a box, cylinder, torus, and annulus, with boundary conditions such that 
various linear operators are separable so that fast direct methods can be applied. The re- 
sulting algorithm is sufficiently efficient that aspect ratios up to F pa 20 can be studied on a 
single-processor workstation over several days. We verified that the code was second-order- 
accurate with respect to the spatial and time resolutions, and that it gave good agreement 
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for the critical Rayleigh number and for the Nusselt number versus Rayleigh number curve 
near onset. 

The most significant feature of our algorithm is the use of a single non-staggered mesh 
for discretizing the equations and fields (velocity, temperature, and pressure). The use of 
a single mesh simplifies the writing and validation of the code, and facilitates adding new 
physical terms like a Coriolis force. The single mesh also allowed the use of fast direct 
methods from the FISHPACK library [24] to solve the Helmholtz and Poisson equations 
associated with the implicit part of each time step. We found that numerical integrations 
of the Boussinesq equations were stable on a single mesh despite results of some previous 
papers that suggested that a non-staggered Navier-Stokes code could be unstable because 
of pressure oscillations. 

Although the algorithm is useful and has been successfully applied to several problems [4, 
5] , there are two ways that the algorithm could be improved for the future study of large- 
aspect-ratio Rayleigh-Benard convection. First is to parallelize the code for a distributed 
memory parallel computer so that aspect ratios comparable to the largest experiments (50 < 
r < 100) could be studied. This is technically straightforward and would involve, first, 
distributing the arrays that represent the fields over the various processors, and, second, 
replacing the fast direct solvers with iterative methods for sparse matrices. Because of 
the simpler data structures and reduced communication overhead associated with a finite- 
difference discretization, the parallelized algorithm will likely be more efficient for simple 
geometries and for simple boundary conditions than the parallel spectral element method 
of Kefs. [13,14]. 

Because time integration algorithms involve sequential steps, parallelizing a code allows 
a larger spatial domain, but not a longer observation time, to be studied for a fixed amount 
of wall-clock time. A second helpful improvement would be to increase the efficiency of 
the time integration method close to the onset of convection so that larger time steps can 
be taken for a given computational effort. A weakness of the operating splitting used in 
most convection codes — finite difference, spectral, and spectral-element — is that the explicit 
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integration of the advection terms imposes a bound on the time step of the form Ce~^/^Ax 
where C is a constant and e = {R — Rc) / Rc is the reduced Rayleigh number. This bound 
is independent of the spatial resolution and diverges less rapidly in the limit e ^ 0"*" than 
the physical time scale, which is proportional to . It would be interesting to explore 
whether a more sophisticated explicit time-stepping technique such as a matrix exponential 
method [40,41] or a fully implicit method [42,43] may succeed in allowing larger time steps 
that are commensurate with the physical time scale while retaining the efficiency of the 
present code. 
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FIGURES 

FIG. 1. (a) Plot of the maximum stable time step as a function of Rayleigh number. The 

Prandtl number is kept constant at o" = 1 (square symbols) and o" = 10 (for cross symbols). The 
cell has aspect ratio Tx = Ty = 2 and periodic sidewalls. The mesh resolution is Ax = 1/16. Small 
random perturbations in the temperature field are used as initial conditions. The simulation is run 
until 20 vertical diffusion times, at which time the Euclidean norm ||T|| is then calculated. The 
value of At such that this norm becomes greater than 10^ is defined as the maximum stable time 
step, (b) Plot of the maximum stable time step as a function of Prandtl number. The Rayleigh 
number was kept constant at R = 2048 (square symbols) and R = 8192 (cross symbols). The same 
aspect ratio, mesh resolution, and initial conditions as in (a) were used. 

FIG. 2. (a) Comparison of numerical (square symbols) and theoretical (cross symbols) 
Nusselt numbers versus Rayleigh number. The numerical values come from time-independent 
two-dimensional nonlinear states (at time t = 12) with periodic sidewalls. The parameters have 
the values = 2, R = 2500, a = 0.71, Ax = Az = 1/16 {N = 16), and At = 0.01. The initial 
state was a small random perturbation of the linear conducting profile. The theoretical values come 
from an asymptotic expansion [34]. (b) Comparison of Nusselt number versus Rayleigh number 
obtained numerically for our algorithm (squares) and for a spectral code of Clever and Busse [35] 
(*). The agreement is better than 3%. 

FIG. 3. (a) Contour lines of the temperature field T{x,z) observed at time t = 12 in a 
two-dimensional box of aspect ratio Tx = 2 with periodic sidewalls. The parameters have val- 
ues R = 2500, a = 0.71, Ax = Az = 1/16 (A^ = 16), and At = 0.01. (b) Contour lines of the 
temperature field observed in a simulation using the same geometry and spatial resolution as in 
(a) but for R = 10^ and time step At = 0.0025. 
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FIG. 4. (a) Time-independent velocity field v(a;, z) at time t = 12 ioi R = 2500, for the same 
geometry and resolution as Fig. 3. The steady state consists of two convection rolls at the critical 
wavenumber Qc = 3.117. (b) Vertical velocity component w{x,z = 0) through the midline of the 
cell, indicating the range of w. 

FIG. 5. (a) Weakly time-dependent temperature contours at the midplane, T{x,y,z = 0), at 
time t = 200 obtained from a three-dimensional box of aspect ratio = Ty = 16 with no-slip 
and insulating boundary conditions, Eqs. (5) and (7). Parameter values are R = 2500, a = 0.71, 
Ax = Ay = Az = 1/8, and At = 0.01. (b): Time-dependent temperature contours at the 
midplane, T{x,y,z = 0), for the same geometry and resolutions as in (a) but for R = 8500 for 
which the rolls arc unstable to the oscillatory instability, which shows up as propagating ripples 
along the rolls. The time-averaged Nusselt number (A'") = 2.27 is larger than for (a), for which the 
value is (AT) = 1.44. 
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TABLES 

TABLE I. Estimated order of convergence ph from Eq. (27) , as a function of the number of mesh 

points N = Tx/h, for a stationary solution of a two-dimensional square box with periodic sidewalls. 
The aspect ratio Tx = 2.016, Rayleigh number R = 1725, and Prandtl number a = 0.71. Results 

■ATC prosculed for for Iho Icrnpcraluro Hold T(x,z) and for the --compoucut of llic volocilv il{x,z). 



N 


Ph for T Ph for w 


16 
32 


1.46 1.43 
1.80 1.79 


TABLE IL Estimated critical Rayleigh number Rc, based on where the growth rate a = (t{R) 
linearly interpolates through zero. The relative error is defined by (Rc — 1708)/1708. 


N 


Rc Relative error 


16 
32 
64 


1693.0 0.9 
1696.6 0.7 
1698.5 0.5 
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Figure 3. 
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Figures 4. 
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Figure 5. 
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